A comparison of JULES-ES-1p0 wave01 members against the original ensemble (wave00).
Wave 01 input parameter sets were picked using History matching to fall within Andy Wiltshires basic constraints on NBP, NPP, cSoil and cVeg stocks at the end of the 20th century. We use 300 of the 500 members, keeping back 2/5ths for emulator validation later.
We answer some basic questions.
What proportion of the new ensemble match AW’s constraints?
How good is a GP emulator? Does it get better overall with the new ensemble members added? In particular, does it get better for those members within the AW constraints?
Does the sensitivity analysis change?
Load libraries, functions and data.
There are no NAs but some relative humidity values are infinite. There are no “low NPP” ensemble members
## [1] 117464.6
## [1] FALSE
## row col
## [1,] 140 9
## [2,] 232 9
## [3,] 249 9
## [4,] 300 9
## [1] Inf Inf Inf Inf
## [1] "rh_lnd_sum"
Global mean for the 20 years at the end of the 20th Century. There is still a significant low bias on cVeg output.
### What proportion of models now fall within Andy’s constraints?
A third! Better than before, but still not great. Pointing at a significant model discrepency in cVeg
Of the 300 members of the wave01 ensemble, 100 pass Andy Wiltshire’s Level 2 constraints.
## [1] 100
Pairs plot of the inputs that pass the constraints with respect to the limits of the original ensemble.
Wave00 in blue and wave01 in red.
Plot Wave00 and Wave01 on top of one another.
We hope that running the new ensemble gives us a better emulator, and allows us to rule out more input space. We particularly hope that the emulator is better for those members that are inside AW’s constraints.
First, we can look at the emulator errors in two cases: The level1a data (a basic carbon cycle), and then with the Wave01 data, which should have similar characteristics. (We should have eliminated really bad simulations, but wave01 is not constrained the data perfectly to be within AW constraints.)
## nbp_lnd_sum npp_nlim_lnd_sum cSoil_lnd_sum cVeg_lnd_sum
## 583 213 317 312
## nbp_lnd_sum npp_nlim_lnd_sum cSoil_lnd_sum cVeg_lnd_sum
## 314 243 243 243
Found the outlier - looks like it’s 440
## integer(0)
The top row shows the leave-one-out prediction accuracy of the original wave00 ensemble, and the lower row the entire wave00 AND wave01 ensemble combined.
We see that the error stats for some of the outputs from wave01 are worse, but there are many more ensemble members that lie within the constraints for wave 01.
“pmae” is “proportional mean absolue error”, which is the mean absolute error expressed as a percentage of the original (minimally constrained) ensemble range in that output.
This gives us an idea of how good the emulator is where it really matters, and as the members are consistent, gives us a fairer idea of whether the emulators have improved with more members.
Good news is, the emulators are more accurate for wave01.
These leave-one-out prediction accuracy plots rank the ensemble members from largest underprediction to largest overprediction using the wave00 predictions. A perfect prediction would appear on the horizontal “zero” line.
Many of the wave01 predictions are closer to the horizontal line, and therefore more accurate predictions.
None of the predictions are outside the uncertainty bounds, which suggests they are overconservative (should be smaller).
Looking at the proportional mean absolute error (pmae), expressed in percent, we can see that it doesn’t improve much for the whole ensemble, but does improve significantly for the subset of ensemble members that fall within AW’s constraints from the first ensemble (marked "_sub").
pmae_wave00 <- lapply(loostats_km_Y_level1a, FUN = function(x) x$pmae )
pmae_wave01 <- lapply(loostats_km_Y_level1a_wave01, FUN = function(x) x$pmae )
pmae_wave00_sub <- lapply(loostats_km_Y_level1a_sub, FUN = function(x) x$pmae )
pmae_wave01_sub <- lapply(loostats_km_Y_level1a_wave01_sub, FUN = function(x) x$pmae )
pmae_table <- cbind(pmae_wave00, pmae_wave01, pmae_wave00_sub, pmae_wave01_sub)
print(pmae_table)
## pmae_wave00 pmae_wave01 pmae_wave00_sub pmae_wave01_sub
## [1,] 5.084658 4.927695 7.138786 4.91346
## [2,] 4.281999 4.007556 4.804586 4.085456
## [3,] 3.597272 3.790159 4.555993 3.834034
## [4,] 4.241615 4.516118 4.81537 3.226458
Calculate the atmospheric growth rate of 1984- 2013 using a simple linear fit
## [1] "correlation agr vs cnbp (all members)"
## [2] "-0.0407011946805412"
## [1] "correlation agr vs cnbp (wave01)" "0.00334447281747124"
##
## Call:
## lm(formula = agr_modern_ens ~ cnbp_modern_ens)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.047413 -0.003802 0.003838 0.006715 0.023090
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.114e-01 5.176e-04 215.173 <2e-16 ***
## cnbp_modern_ens -1.379e-09 1.519e-09 -0.908 0.364
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01155 on 497 degrees of freedom
## Multiple R-squared: 0.001657, Adjusted R-squared: -0.0003522
## F-statistic: 0.8247 on 1 and 497 DF, p-value: 0.3643
##
## Call:
## lm(formula = agr_modern_wave01 ~ cnbp_modern_wave01)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4639 -4639 -4639 -4639 1382335
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.639e+03 4.639e+03 1.000 0.318
## cnbp_modern_wave01 5.475e-29 9.483e-28 0.058 0.954
##
## Residual standard error: 80210 on 298 degrees of freedom
## Multiple R-squared: 1.119e-05, Adjusted R-squared: -0.003344
## F-statistic: 0.003333 on 1 and 298 DF, p-value: 0.954
Interannual variability and cumulative NBP
(correlations are close to zero, especially in the later wave)
## [1] 0.03223847
## [1] 0.003344484
Using Atmospheric Growth Rate as an example, how close can we get the model to observations? Can we do better than standard? What are the trade offs of doing so? How does getting close in AGR affect performance in other outputs?
We’ve established that most of the original ensemble have an ME/MAE/RMSE larger than the standard run. More (but few) of the wave01 perform better than standard.
A map of the 2D projections of parameter space where the ensemble member performs better than standard.
The blue part is the first wave, and not subject to constraint so may be removed in the second wave (wave01).
Having trouble fitting RMSE, to trying mean error.
Why is there an odd collection at just under 1?
## [1] 100
This next pairs plot looks at all the ensemble members that have a better mean atmospheric growth error than standard.
This next plot looks at all the ensemble members that have a better mean atmospheric growth error than standard AND pass the level 2 constraints.
The number is small (41/300), but the ensemble members seem spread across parameter space.
##
## optimisation start
## ------------------
## * estimation method : MLE
## * optimisation method : BFGS
## * analytical gradient : used
## * trend model : ~alpha_io + a_wl_io + bio_hum_cn + b_wl_io + dcatch_dlai_io +
## dqcrit_io + dz0v_dh_io + f0_io + fd_io + g_area_io + g_root_io +
## g_wood_io + gs_nvg_io + hw_sw_io + kaps_roth + knl_io + lai_max_io +
## lai_min_io + lma_io + n_inorg_turnover + nmass_io + nr_io +
## retran_l_io + retran_r_io + r_grow_io + rootd_ft_io + sigl_io +
## sorp + tleaf_of_io + tlow_io + tupp_io + l_vg_soil
## * covariance model :
## - type : matern5_2
## - nugget : NO
## - parameters lower bounds : 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
## - parameters upper bounds : 1.978915 1.984434 1.985775 1.678886 1.961948 1.960223 1.976754 1.774573 1.971703 1.983945 1.989915 1.994414 1.970738 1.986029 1.924591 1.992288 1.984011 1.960324 1.994092 1.983781 1.953327 1.972882 1.972018 1.982383 1.986984 1.982191 1.988709 1.9928 1.987255 1.995015 1.96482 1.952566
## - best initial criterion value(s) : -93.32195
##
## N = 32, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate 0 f= 93.322 |proj g|= 1.3048
## At iterate 1 f = 92.672 |proj g|= 1.9402
## At iterate 2 f = 92.632 |proj g|= 1.903
## At iterate 3 f = 92.573 |proj g|= 1.9385
## At iterate 4 f = 92.455 |proj g|= 1.9421
## At iterate 5 f = 92.236 |proj g|= 1.941
## At iterate 6 f = 90.691 |proj g|= 1.6648
## At iterate 7 f = 90.435 |proj g|= 1.9429
## At iterate 8 f = 90.217 |proj g|= 1.6173
## At iterate 9 f = 90.103 |proj g|= 1.5773
## At iterate 10 f = 89.956 |proj g|= 1.5013
## At iterate 11 f = 89.75 |proj g|= 1.9419
## At iterate 12 f = 89.443 |proj g|= 1.944
## At iterate 13 f = 89.373 |proj g|= 1.9435
## At iterate 14 f = 89.257 |proj g|= 1.8349
## At iterate 15 f = 89.23 |proj g|= 1.9439
## At iterate 16 f = 89.164 |proj g|= 1.826
## At iterate 17 f = 88.923 |proj g|= 1.8345
## At iterate 18 f = 88.43 |proj g|= 1.8078
## At iterate 19 f = 88.143 |proj g|= 1.7868
## At iterate 20 f = 87.359 |proj g|= 1.9417
## At iterate 21 f = 87.294 |proj g|= 1.6155
## At iterate 22 f = 87.241 |proj g|= 1.6026
## At iterate 23 f = 87.214 |proj g|= 1.9408
## At iterate 24 f = 87.19 |proj g|= 1.9414
## At iterate 25 f = 87.161 |proj g|= 1.9416
## At iterate 26 f = 87.128 |proj g|= 1.9416
## At iterate 27 f = 86.934 |proj g|= 1.945
## At iterate 28 f = 86.756 |proj g|= 1.9439
## At iterate 29 f = 86.51 |proj g|= 1.4558
## At iterate 30 f = 86.455 |proj g|= 1.7678
## At iterate 31 f = 86.334 |proj g|= 0.77311
## At iterate 32 f = 86.311 |proj g|= 1.9419
## At iterate 33 f = 86.273 |proj g|= 1.942
## At iterate 34 f = 85.9 |proj g|= 1.9424
## At iterate 35 f = 85.798 |proj g|= 1.6463
## At iterate 36 f = 85.786 |proj g|= 1.6456
## At iterate 37 f = 85.753 |proj g|= 1.6368
## At iterate 38 f = 85.71 |proj g|= 1.9422
## At iterate 39 f = 85.665 |proj g|= 1.7838
## At iterate 40 f = 85.662 |proj g|= 1.7837
## At iterate 41 f = 85.639 |proj g|= 0.95615
## At iterate 42 f = 85.621 |proj g|= 0.65895
## At iterate 43 f = 85.486 |proj g|= 1.4736
## At iterate 44 f = 85.376 |proj g|= 1.9418
## At iterate 45 f = 85.121 |proj g|= 1.9417
## At iterate 46 f = 85.116 |proj g|= 0.47806
## At iterate 47 f = 84.971 |proj g|= 0.66045
## At iterate 48 f = 84.962 |proj g|= 1.9416
## At iterate 49 f = 84.934 |proj g|= 1.942
## At iterate 50 f = 84.918 |proj g|= 1.9418
## At iterate 51 f = 84.896 |proj g|= 0.78518
## At iterate 52 f = 84.886 |proj g|= 0.67606
## At iterate 53 f = 84.881 |proj g|= 0.6382
## At iterate 54 f = 84.874 |proj g|= 0.3609
## At iterate 55 f = 84.872 |proj g|= 1.6025
## At iterate 56 f = 84.866 |proj g|= 1.9411
## At iterate 57 f = 84.848 |proj g|= 1.7796
## At iterate 58 f = 84.827 |proj g|= 0.98429
## At iterate 59 f = 84.821 |proj g|= 0.35486
## At iterate 60 f = 84.82 |proj g|= 0.3659
## At iterate 61 f = 84.818 |proj g|= 1.3651
## At iterate 62 f = 84.813 |proj g|= 1.9415
## At iterate 63 f = 84.802 |proj g|= 1.9419
## At iterate 64 f = 84.776 |proj g|= 1.9423
## At iterate 65 f = 84.762 |proj g|= 1.9421
## At iterate 66 f = 84.736 |proj g|= 0.47123
## At iterate 67 f = 84.733 |proj g|= 0.8154
## At iterate 68 f = 84.732 |proj g|= 0.471
## At iterate 69 f = 84.73 |proj g|= 0.47152
## At iterate 70 f = 84.724 |proj g|= 0.47237
## At iterate 71 f = 84.704 |proj g|= 0.47447
## At iterate 72 f = 84.663 |proj g|= 0.47812
## At iterate 73 f = 84.618 |proj g|= 0.47998
## At iterate 74 f = 84.585 |proj g|= 1.9416
## At iterate 75 f = 84.572 |proj g|= 1.9414
## At iterate 76 f = 84.554 |proj g|= 0.47151
## At iterate 77 f = 84.462 |proj g|= 1.1364
## At iterate 78 f = 84.454 |proj g|= 1.2454
## At iterate 79 f = 84.451 |proj g|= 1.3069
## At iterate 80 f = 84.45 |proj g|= 1.3387
## At iterate 81 f = 84.45 |proj g|= 1.6226
## At iterate 82 f = 84.449 |proj g|= 1.9409
## At iterate 83 f = 84.445 |proj g|= 1.2623
## At iterate 84 f = 84.435 |proj g|= 0.85365
## At iterate 85 f = 84.422 |proj g|= 1.6125
## At iterate 86 f = 84.42 |proj g|= 0.33685
## At iterate 87 f = 84.415 |proj g|= 0.32882
## At iterate 88 f = 84.411 |proj g|= 0.31766
## At iterate 89 f = 84.405 |proj g|= 0.36348
## At iterate 90 f = 84.383 |proj g|= 0.3772
## At iterate 91 f = 84.361 |proj g|= 0.39268
## At iterate 92 f = 84.257 |proj g|= 1.4756
## At iterate 93 f = 84.225 |proj g|= 1.9415
## At iterate 94 f = 84.207 |proj g|= 1.941
## At iterate 95 f = 84.203 |proj g|= 1.941
## At iterate 96 f = 84.197 |proj g|= 1.941
## At iterate 97 f = 84.185 |proj g|= 0.46746
## At iterate 98 f = 84.181 |proj g|= 0.45837
## At iterate 99 f = 84.164 |proj g|= 1.9422
## At iterate 100 f = 84.149 |proj g|= 1.9414
## At iterate 101 f = 84.144 |proj g|= 0.49707
## At iterate 102 f = 84.144 |proj g|= 0.49703
## At iterate 103 f = 84.14 |proj g|= 0.49651
## At iterate 104 f = 84.134 |proj g|= 0.67445
## At iterate 105 f = 84.13 |proj g|= 0.54398
## At iterate 106 f = 84.106 |proj g|= 1.3705
## At iterate 107 f = 84.071 |proj g|= 1.6366
## At iterate 108 f = 84.052 |proj g|= 1.9414
## At iterate 109 f = 84.036 |proj g|= 1.9416
## At iterate 110 f = 84.017 |proj g|= 1.9412
## At iterate 111 f = 84.011 |proj g|= 1.852
## At iterate 112 f = 84.01 |proj g|= 0.69667
## At iterate 113 f = 84.009 |proj g|= 0.70693
## At iterate 114 f = 84.006 |proj g|= 0.73004
## At iterate 115 f = 83.994 |proj g|= 0.78937
## At iterate 116 f = 83.992 |proj g|= 0.80705
## At iterate 117 f = 83.967 |proj g|= 0.90411
## At iterate 118 f = 83.926 |proj g|= 1.7211
## At iterate 119 f = 83.908 |proj g|= 1.9416
## At iterate 120 f = 83.868 |proj g|= 1.9415
## At iterate 121 f = 83.848 |proj g|= 1.9416
## At iterate 122 f = 83.825 |proj g|= 0.78089
## At iterate 123 f = 83.786 |proj g|= 0.73467
## At iterate 124 f = 83.741 |proj g|= 0.70987
## At iterate 125 f = 83.734 |proj g|= 1.9411
## At iterate 126 f = 83.722 |proj g|= 1.9411
## At iterate 127 f = 83.653 |proj g|= 1.941
## At iterate 128 f = 83.556 |proj g|= 1.9411
## At iterate 129 f = 83.525 |proj g|= 1.9413
## At iterate 130 f = 83.524 |proj g|= 1.5405
## At iterate 131 f = 83.48 |proj g|= 1.4828
## At iterate 132 f = 83.449 |proj g|= 0.73528
## At iterate 133 f = 83.432 |proj g|= 0.67389
## At iterate 134 f = 83.425 |proj g|= 0.23426
## At iterate 135 f = 83.424 |proj g|= 0.14269
## At iterate 136 f = 83.424 |proj g|= 0.14252
## At iterate 137 f = 83.424 |proj g|= 0.14272
## At iterate 138 f = 83.423 |proj g|= 0.14133
## At iterate 139 f = 83.415 |proj g|= 1.942
## At iterate 140 f = 83.41 |proj g|= 1.7477
## At iterate 141 f = 83.407 |proj g|= 1.1435
## At iterate 142 f = 83.406 |proj g|= 0.21509
## At iterate 143 f = 83.406 |proj g|= 0.15963
## At iterate 144 f = 83.404 |proj g|= 0.15344
## At iterate 145 f = 83.402 |proj g|= 0.25745
## At iterate 146 f = 83.398 |proj g|= 0.25647
## At iterate 147 f = 83.39 |proj g|= 0.34118
## At iterate 148 f = 83.39 |proj g|= 0.33371
## At iterate 149 f = 83.39 |proj g|= 0.34581
## At iterate 150 f = 83.39 |proj g|= 0.84181
## At iterate 151 f = 83.39 |proj g|= 0.42181
## At iterate 152 f = 83.39 |proj g|= 0.35403
## At iterate 153 f = 83.389 |proj g|= 0.34286
## At iterate 154 f = 83.389 |proj g|= 0.36132
## At iterate 155 f = 83.387 |proj g|= 0.47538
## At iterate 156 f = 83.384 |proj g|= 0.3362
## At iterate 157 f = 83.382 |proj g|= 0.35209
## At iterate 158 f = 83.38 |proj g|= 0.42829
## At iterate 159 f = 83.38 |proj g|= 0.3877
## At iterate 160 f = 83.38 |proj g|= 0.38009
## At iterate 161 f = 83.379 |proj g|= 0.34778
## At iterate 162 f = 83.378 |proj g|= 0.25188
## At iterate 163 f = 83.378 |proj g|= 0.50895
## At iterate 164 f = 83.377 |proj g|= 0.17723
## At iterate 165 f = 83.376 |proj g|= 0.11848
## At iterate 166 f = 83.376 |proj g|= 0.078262
## At iterate 167 f = 83.376 |proj g|= 0.12043
## At iterate 168 f = 83.376 |proj g|= 0.061317
## At iterate 169 f = 83.376 |proj g|= 1.2626
## At iterate 170 f = 83.376 |proj g|= 0.04529
## At iterate 171 f = 83.376 |proj g|= 0.029058
## At iterate 172 f = 83.376 |proj g|= 0.029054
##
## iterations 172
## function evaluations 198
## segments explored during Cauchy searches 175
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 26
## norm of the final projected gradient 0.029054
## final function value 83.3757
##
## F = 83.3757
## final value 83.375724
## converged
## [1] 21.181
This pairs plot shows the 2d and marginal density of emulated input points where the emulated atmospheric growth is closer to the observations than the standard model.
This technique might provide a useful set of points for optimising the model (at least to atmospheric growth).
Next, check emulators of all the other outputs and apply the constraints to them. See how the constraints change.
## [1] 12.947
Emulated members passing level2 constraints AND having lower error in atmospheric growth than standard.
Red point indicates the standard input.
## [1] 5.402